#!/usr/bin/perl -wT
use strict;
use CGI::Carp qw(fatalsToBrowser warningsToBrowser);
use DBI;
use DBD::Pg qw(:pg_types);
use Data::Dumper;
use GD;
use GD::Text;
use POSIX qw(ceil floor);

# Open a database connection:
my $dbh = DBI->connect("dbi:Pg:dbname=protein", "username", "password",{ShowErrorStatement => 1,PrintError => 1});
if (!defined $dbh) {die($DBI::errstr);}

# Database variables:
my $rc;
my $sth;

sub dbfuncs {
  my $vars;
  $vars->{dbfuncs}=1;
  $vars->{getstudylist}=\&getstudylist;
  $vars->{getstudy}=\&getstudy;
  $vars->{getexperimentlist}=\&getexperimentlist;
  $vars->{getexperiment}=\&getexperiment;
  $vars->{getspectra_peptidelist}=\&getspectra_peptidelist;
  $vars->{getspectra_peptide}=\&getspectra_peptide;
  $vars->{getprotein_quantificationlist}=\&getprotein_quantificationlist;
  $vars->{getprotein_quantification}=\&getprotein_quantification;
  $vars->{getproteinlist}=\&getproteinlist;
  $vars->{getprotein}=\&getprotein;
  $vars->{getacq_operatorlist}=\&getacq_operatorlist;
  $vars->{getacq_operator}=\&getacq_operator;
  $vars->{getacq_mass_speclist}=\&getacq_mass_speclist;
  $vars->{getacq_mass_spec}=\&getacq_mass_spec;
  $vars->{getacq_analytical_columnlist}=\&getacq_analytical_columnlist;
  $vars->{getacq_analytical_column}=\&getacq_analytical_column;
  $vars->{getacq_mobile_phase_alist}=\&getacq_mobile_phase_alist;
  $vars->{getacq_mobile_phase_a}=\&getacq_mobile_phase_a;
  $vars->{getacq_mobile_phase_blist}=\&getacq_mobile_phase_blist;
  $vars->{getacq_mobile_phase_b}=\&getacq_mobile_phase_b;
  $vars->{getsearch_enginelist}=\&getsearch_enginelist;
  $vars->{getsearch_engine}=\&getsearch_engine;
  $vars->{getsearch_proteinbanklist}=\&getsearch_proteinbanklist;
  $vars->{getsearch_proteinbank}=\&getsearch_proteinbank;
  $vars->{wrapprotein}=\&wrapprotein;
  return $vars;
}

sub wrapprotein {
  my ($study,$experiment,$protein_quantification,$spectra_peptide)=@_;
  #return "<tt>study=$study,experiment=$experiment,protein_quantification=$protein_quantification,spectra_peptide=$spectra_peptide</tt>";

  # Find the protein sequence for this pq:
  my $r=$dbh->selectall_hashref("select pq.id,p.sequence from protein as p, protein_quantification as pq where pq.id=? and p.id=pq.protein_id","id",undef,$protein_quantification);
  if (!defined $r) {return "The protein for this protein_quantification was not found";}
  my $sequence=$r->{$protein_quantification}{sequence};

  # Find the spectra/peptides that go with this pq:
  $r=$dbh->selectall_hashref("select id,peptide_sequence from spectra_peptide
  where protein_quantification_id=?","id",undef,$protein_quantification);
  if (!defined $r) {return "That protein had no supporting spectra/peptides";}
  my %others=map {$_ => $r->{$_}{peptide_sequence}} keys %$r;

  # Find the current peptide:
  my $current="";
  if (defined $spectra_peptide and $spectra_peptide ne "") {
    $r=$dbh->selectall_hashref("select id,peptide_sequence from spectra_peptide
  where id=?","id",undef,$spectra_peptide);
    $current=$r->{$spectra_peptide}{peptide_sequence};
  }

  return fixit($sequence,$current,%others);
}

#############################################################################
# Wrap a long sequence by inserting spaces every 'n' characters:
sub wrapit {
  my ($string)=$_[0];
  if (!defined $string) {return undef;}
  $string=~s/(.{80})/$1<BR\/>\n/gs;
  return $string;
}

#############################################################################
# Wrap a long sequence and highlight a subsequence:
# String: protein sequence
# Substring: peptide sequence
# Others: key=peptide id, value=peptide sequence
sub fixit {
  my ($string,$substring,%others)=@_;
  my @link;
  my $m;
  my $n;
  my $baselink=$ENV{'REQUEST_URI'};
  #return "<tt>".Dumper(\%ENV)."</tt>";
  #return "<tt>$baselink</tt>";
  $baselink=~s,spectra_peptide-\d+/,,g;
  my $color="";  # Color: blank or 'allmatch'
  my $alink="";  # A-link: link of previous letter
  my $blink="";  # B-link: link of this letter
  my $newstring="";  # The new constructed string
  my $letter="";  # The letter

  $substring=uc($substring);
  $string=uc($string);  # Uppercase both strings

  # Calculate links for each letter for the '%others' list:
  foreach (keys %others) {  # For each peptide,
    if (($m=index($string,$others{$_}))>=0) {  # If the peptide is found in this protein,
      for ($n=$m;$n<$m+length($others{$_});$n++) {  # Attach a link to each letter
        $link[$n]="${baselink}spectra_peptide-$_/";  # Last link wins
      }
    }
  }

  # Highlight the current peptide:
  if (defined($substring) and $substring ne "") {
    if (($m=index($string,$substring))>=0) {  # If the peptide is found in this protein,
      for ($n=$m;$n<$m+length($substring);$n++) {  # Attach a link to each letter
	$link[$n]=".";                 # Last link wins
      }
    }
  }


  # Construct new string:
  for ($n=0;$n<length($string);$n++) {  # For each letter of original string
    $blink=$link[$n];
    if (!defined $blink) {$blink="";}
    $letter=substr($string,$n,1); # Get the current letter
    if ($n%60==0 and $n!=0) { $newstring.="<BR/>\n"; }
    if ($blink ne $alink) {       # If link is different
      if ($alink ne "") { $newstring.="</span></a>";}  # End previous link
      if ($blink ne "") {
        $newstring.="<a href='$blink'><span class=";
	$newstring.=($blink=~m/[.]/)?"currentmatch":"allmatch";
	$newstring.=">";
      } # Start of new link
      $alink=$blink;              # this link becomes the new 'previous' link
    }
    $newstring.=$letter;   # Write the current letter.
  }
  if ($alink ne "") {$newstring.="</a>";}
  return $newstring;
}


sub getstudylist {return getit("study","where 1=1 or 1=?",$_[0]);}
sub getexperimentlist {return getit("experiment","where study_id=?",$_[0]);}

sub getspectra_peptidelist {return $dbh->selectcol_arrayref("select s.id from spectra_peptide as s where s.experiment_id=? order by s.peptide_sequence",undef,$_[0]);}

sub getprotein_quantificationlist {return $dbh->selectcol_arrayref("select pq.id from protein_quantification as pq, protein as p where experiment_id=? and pq.protein_id = p.id order by p.accession_code",undef,$_[0]);}

sub getproteinlist {return getit("protein","where 1=1 or 1=?",$_[0]);}
sub getacq_operatorlist {return getit("acq_operator","where 1=1 or ''=?","");}
sub getacq_mass_speclist {return getit("acq_mass_spec","where 1=1 or ''=?","");}
sub getacq_analytical_columnlist {return getit("acq_analytical_column","where 1=1 or ''=?","");}
sub getacq_mobile_phase_alist {return getit("acq_mobile_phase_a","where 1=1 or ''=?","");}
sub getacq_mobile_phase_blist {return getit("acq_mobile_phase_b","where 1=1 or ''=?","");}
sub getsearch_enginelist {return getit("search_engine","where 1=1 or ''=?","");}
sub getsearch_proteinbanklist {return getit("search_proteinbank","where 1=1 or ''=?","");}

sub getstudy {return getit("study","",$_[0]);}
sub getexperiment {return getit("experiment","",$_[0]);}
sub getspectra_peptide {return getit("spectra_peptide","",$_[0]);}
sub getprotein_quantification {return getit("protein_quantification","",$_[0]);}
sub getprotein {return getit("protein","",$_[0]);}
sub getacq_operator {return getit("acq_operator","",$_[0]);}
sub getacq_mass_spec {return getit("acq_mass_spec","",$_[0]);}
sub getacq_analytical_column {return getit("acq_analytical_column","",$_[0]);}
sub getacq_mobile_phase_a {return getit("acq_mobile_phase_a","",$_[0]);}
sub getacq_mobile_phase_b {return getit("acq_mobile_phase_b","",$_[0]);}
sub getsearch_engine {return getit("search_engine","",$_[0]);}
sub getsearch_proteinbank {return getit("search_proteinbank","",$_[0]);}

sub getit {
  my ($what,$type,$arg) = @_;
  my $sql;
  my @stuff;

  while (ref($arg) eq "REF") {$arg=$$arg;}
  if (ref($arg) eq "ARRAY") {@stuff=@$arg;}
  if (ref($arg) eq "") {@stuff=($arg);}

  if ($type ne "") {
    # Get a list of ids:
    # select id from foo where 1 = 1 or 1 = ?
    # select id from bar where experiment_id = ?
    # select id from experiment where study_id = ?
    $sql="select id from $what $type order by id";
    return $dbh->selectcol_arrayref($sql,undef,@stuff);
  }

  # Get a row of data:
  # select * from foo where id = ?
  if ($stuff[0] eq "") {return undef;}
  $sql="select * from $what where id=?";
  return $dbh->selectrow_hashref($sql,undef,@stuff);
}

sub getoptions { 
  my $r=$dbh->selectall_hashref("select * from ".($_[0])." order by id","id");
  if (!defined $r) {
    die("getoptions called with table $_[0], but got no result.");
  }
  my %result=%{$r};
  my %r2=map {$_ => $result{$_}{name}} keys %result;
  return \%r2;
}

# Do a column query without the caller having to connect to the db.
sub getlist {
  #warn("getlist called: sql=$_[0]");
  my $r=$dbh->selectcol_arrayref($_[0]);
  #warn(Dumper($r)." ");
  return $r;
}

# Do a query without the caller having to connect to the db.
sub selectall_arrayref {
#  warn("getarray called: sql=$_[0]");
  my $r=$dbh->selectall_arrayref(shift @_,undef,@_);
#  warn(Dumper($r)." ");
#  exit 0;
  return $r;
}

# Graphics:

# Read in data:
#my (@x,@y,@ann);
#open(IN,"<s1.txt") or die("s1.txt: $!");
#while (<IN>) {
#  chomp;
#  if (m/^([\d.]+) ([\d.]+)/) {
#    push @x,$1;
#    push @y,$2;
#    push @ann,"";
#  }
#}
#close(IN);

## Label the high points.
#my $n;
#my $maxy=$y[0];
#for ($n=1;$n<@y;$n++) { if ($y[$n]>$maxy) {$maxy=$y[$n];} }
#for ($n=1;$n<@y;$n++) { if ($y[$n]>0.5*$maxy) {$ann[$n]=$x[$n];} }

#binmode STDOUT;
#print graph(\@x,\@y,\@ann,500,400);

# Amino acid masses:
my %aamass=(
  "A"=> 71.03712,
  "R"=> 156.10112,
  "N"=> 114.04293,
  "D"=> 115.02695,
  "C"=> 103.00919,
  "E"=> 129.04264,
  "Q"=> 128.05858,
  "G"=> 57.02147,
  "H"=> 137.05891,
  "I"=> 113.08407,
  "L"=> 113.08407,
  "K"=> 128.09497,
  "M"=> 131.04049,
  "F"=> 147.06842,
  "P"=> 97.05277,
  "S"=> 87.03203,
  "T"=> 101.04768,
  "W"=> 186.07932,
  "Y"=> 163.06333,
  "V"=> 99.06842,
);

sub simulatespectra {
  my $in=$_[0];
  my @masses;
  my @intensities;
  if (!defined $in or $in eq "") {die("simulatespectra called with no parameters");}
  my $m=0;
  my $i=0;
  my $r="";
  for (my $x=0;$x<length($in);$x++) {
    $m+=$aamass{substr($in,$x,1)};
    $i=rand(80)+30;
    $r.="$m $i\n";
  }
  return $r;
}

# Returns a PNG graph of the spectra
sub graphspectra {
  my @x=@{$_[0]}; # Masses
  my @y=@{$_[1]}; # Intensities
  my @ann=@{$_[2]}; # Annotations
  # How wide and high should this graph be?
  my $xwidth=(defined($_[3]))?$_[3]:500;
  my $yheight=(defined($_[4]))?$_[4]:500;
  my $simulation=$_[5]; # Is this a simulation?

  my ($m,$n,$o);
  my $high=textheight("1243");

  # Set up the image, draw the box
  my $im=new GD::Image($xwidth,$yheight);
  my $white=$im->colorAllocate(255,255,255);
  my $black=$im->colorAllocate(0,0,0);
  my $blue=$im->colorAllocate(0,0,192);
  my $grey=$im->colorAllocate(192,192,192);
  $im->fill(0,0,$white);
  $im->rectangle(0.1*$xwidth,0.95*$yheight,0.95*$xwidth,0.05*$yheight,$blue);
  if (defined $simulation && $simulation == 1) {
    for ($m=0;$m<1;$m+=0.2) {
      for ($n=0;$n<1;$n+=0.2) {
	$im->string(gdGiantFont,($m+rand(0.1))*$xwidth,($n+rand(0.2))*$yheight,"Simulated",$grey);
      }
    }
  }

  # Pick arbitraty min and max values to start
  my ($minx,$miny,$maxx,$maxy)=($x[0],$y[0],$x[0],$y[0]);

  # Error checking
  if (!defined $minx) {die("Nothing was read?");}

  # Find min + max
  foreach (@x) {if($_>$maxx){$maxx=$_;};if($_<$minx){$minx=$_;}}
  foreach (@y) {if($_>$maxy){$maxy=$_;};if($_<$miny){$miny=$_;}}

  # Give graph a little headroom (20%)
  $maxy*=1.2;

  # If minimum y value isn't close to 0 or negative, the graph will look bad
  if ($miny>0) {$miny=0;}

  # Find range of x and y values
  my $rangex=$maxx-$minx;
  my $rangey=$maxy-$miny;

  # Round off the scale in X
  my $scalex=$rangex;
  $scalex=~s/[\d]/0/g;
  $scalex=~s/^0/1/;
  if ($rangex=~m/^1/) {$scalex/=5;}

  # Round off the scale in Y
  my $scaley=$rangey;
  $scaley=~s/\d/0/g;
  $scaley=~s/^0/1/;
  if ($rangey=~m/^1/) {$scaley/=5;}

  # Find the start and end values
  my ($startx,$endx);
  my ($starty,$endy);
  $startx=floor($minx/$scalex)*$scalex;
  $starty=floor($miny/$scaley)*$scaley;
  $endx=ceil($maxx/$scalex)*$scalex;
  $endy=ceil($maxy/$scaley)*$scaley;

  # Put tickmarks on the x axis:
  for ($n=$startx;$n<=$endx;$n+=$scalex) {
    $m=scale($n,$startx,$endx,0.1*$xwidth,0.95*$xwidth);
    $im->line($m,0.95*$yheight,$m,0.95*$yheight+4,$blue);
    $im->string(gdSmallFont,$m-textwidth($n)/2,0.95*$yheight+4,$n,$blue);
  }
  # Put tickmarks on the y axis:
  for ($n=$starty;$n<=$endy;$n+=$scaley) {
    $m=scale($n,$starty,$endy,0.95*$yheight,0.05*$yheight);
    $im->line(0.1*$xwidth-4,$m,0.1*$xwidth,$m,$blue);
    $im->string(gdSmallFont,0.1*$xwidth-4-textwidth($n),$m-$high/2,$n,$blue);
  }

  # Putting the m/z and Intensity labels on the axes:
  $im->string(gdSmallFont,0.75*$xwidth,$yheight-$high,"m/z",$blue);
  $im->stringUp(gdSmallFont,0,0.5*$yheight,"Intensity",$blue);

  # Draw the data:
  for ($n=0;$n<@x;$n++) {
    $m=scale($x[$n],$startx,$endx,0.1*$xwidth,0.95*$xwidth);
    $o=scale($y[$n],$starty,$endy,0.95*$yheight,0.05*$yheight);
    $im->line($m,0.95*$yheight,$m,$o,$black);
    if (defined $ann[$n] and $ann[$n] ne "") {
      $im->line($m,$o,$m-4,$o-2,$black);
      $im->stringUp(gdSmallFont,$m-$high,$o-4,$ann[$n],$black);
    }
  }

  # Return the image:
  return $im->png;
}

sub textwidth { # How wide is string 'xyz'?
  my $in=$_[0];
  my $gd_text=GD::Text->new(text=>$in,font=>gdSmallFont);
  return $gd_text->get('width');
}

sub textheight { # How high is string 'xyz'?
  my $in=$_[0];
  my $gd_text=GD::Text->new(text=>$in,font=>gdSmallFont);
  return $gd_text->get('height');
}

sub scale { # Convert a value to new coords.
  my ($in,$minin,$maxin,$minout,$maxout)=@_;
  my $percent=($in-$minin)/($maxin-$minin);
  return $percent*($maxout-$minout)+$minout;
}

1;
